
%%%%%%%%%%%%%%%%%%%
%%% Calibration %%%
%%%%%%%%%%%%%%%%%%%

clc
clear all
tic

% variable set-up

params.alpha=0.5;
params.eta=6;
params.r=0.02;
params.pi=0.001;

guess_cal=[0.024 0.429 0.114 0.091 17.431 454.703 164.026 839.467];

params.weight=[1 1 1 1 1 1 1 1 1];
params.weight=params.weight./sum(params.weight);

iter=0;
dist=1;
while dist>0.001 % average percent deviation 0.1%
    iter=iter+1
% local min
    [xx_cal,fval_cal] = fminsearchbnd(@(x) calibration_function(x,params),guess_cal,[0 0 0 0 0 0 0 0],[1 1 1 1 10^4 10^4 10^4 10^4]);
    dist=0;
    for i=1:8
        dist=dist+(1/8)*(abs(xx_cal(i)-guess_cal(i))/abs(guess_cal(i)));
    end
    dist
    guess_cal=xx_cal;
end

% momdents calculation

    params.lambdaE=xx_cal(1);
    params.lambdaS=xx_cal(2);
    params.lambdaI=xx_cal(3);
    params.lambdaC=xx_cal(4);
    params.f=xx_cal(5);
    params.xiE=xx_cal(6);
    params.xiI=xx_cal(7);
    params.lambdae=params.lambdaE;
    params.xie=xx_cal(8);

    guess_sim=0.1*ones(10,1);
    options = optimset('MaxFunEvals',10^5,'MaxIter',10^5,'TolFun',1e-30);
    [xx_sim,fval_sim] = fminsearchbnd(@(x) equation_difference(x,params),guess_sim,[0 0 0 0 0 0 0 0 0 0],[],options);
    
    alpha=params.alpha;
    pi=params.pi;
    eta=params.eta;
    r=params.r;
    lambdaE=params.lambdaE;
    lambdaS=params.lambdaS;
    lambdaI=params.lambdaI;
    lambdaC=params.lambdaC;
    lambdae=params.lambdae;
    f=params.f;
    xiE=params.xiE;
    xiI=params.xiI;
    xie=params.xie;
    
    xI=xx_sim(2);
    xE=xx_sim(5);
    xe=xx_sim(8);
    tau=xE+xe;
    theta=xe/(xe+xE+xI);
    B=xx_sim(6);
    gq=xx_sim(7);
    RD_all=xiI*eta^(-1)*(xI)^(1/(1-alpha))+xiE*eta^(-1)*(xE)^(1/(1-alpha))+xie*eta^(-1)*(xe)^(1/(1-alpha))+B*xI*eta^(-1)*pi^(-1)*(1-f);
    RD_new=xie*eta^(-1)*(xe)^(1/(1-alpha));
    
    moment1=xI*(eta-1)*log(1-lambdaC)+(xE+xe)*(eta-1)*log(1-lambdaS); % target -0.1506
    moment2=((xE+xe)*log(1-lambdaS))/(xI*log(1-lambdaC)); % target 1.4858
    moment3=xE*((1+lambdaE)^(eta-1)+(1-lambdaS)^(eta-1)-2)+xI*((1+lambdaI)^(eta-1)+(1-lambdaC)^(eta-1)-2)+xe*((1+lambdae)^(eta-1)+(1-lambdaS)^(eta-1)-2); % target 0.016
    moment4=xe*(1+lambdae)^(eta-1); % target 0.0114
    moment5=xI/xE; % target 6
    moment6=RD_all; % target 0.0988
    moment7=xI+xE; % target 0.145
    moment8=xe; % target 0.0108
    moment9=RD_new; % target 0.0115
    
    moments=[moment1 moment2 moment3 moment4 moment5 moment6 moment7 moment8 moment9]
    
    weight=params.weight;
    yyy=weight(1)*(abs(moment1+0.1506)/0.1506) + weight(2)*(abs(moment2-1.4858)/1.4858) + weight(3)*(abs(moment3-0.016)/0.016) + weight(4)*(abs(moment4-0.0114)/0.0114) + weight(6)*(abs(moment6-0.0988)/0.0988) + weight(7)*(abs(moment7-0.145)/0.145) + weight(8)*(abs(moment8-0.0108)/0.0108) + weight(9)*(abs(moment9-0.0115)/0.0115);

    fval_cal
    fval_sim
    
    xx_eqm=[xI xE xe gq]
   
    save calibration_outcome.mat xx_cal xx_eqm moments
    iter
    toc